#include <stdio.h>
#include <string.h>
#include "io_charmm_inp.h"
#include <ctype.h>
#include "esmd_types.h"
#include "list.h"
#include "memory_cpp.hpp"
#include "list_cpp.hpp"
//One should follow https://www.charmm.org/charmm/documentation/by-version/c42b1/params/doc/parmfile/ to write a charmm inp file parser. DO NOT TRUST GOOGLE's RESULT
enum inp_sections {
  PREAMBLE,
  BOND,
  ANGLE,
  TORI,
  IMPR,
  CMAP,
  NONB,
  HBOND
};
#define INP_MAX_LINE_LEN 65536
INLINE char *inp_trim(char *s) {
  char *ret = s;
  while (isspace(*ret) && *ret)
    ret++;
  if (*ret == '!')
    *ret = 0;
  char *end = ret;
  while (*end && *end != '!')
    end++;
  while (end > ret && isspace(*(end - 1)))
    end--;
  if (end > ret)
    *end = 0;
  return ret;
}

charmm_rawparam_t *read_inp(FILE *finp) {
  charmm_rawparam_t *ffpar = esmd::allocate<charmm_rawparam_t>(1, "charmm parameter");
  esmd::list<charmm_raw_bond_param_t> bond_list("charmm bond parameters");
  esmd::list<charmm_raw_angle_param_t> angle_list("charmm angle parameters");
  esmd::list<charmm_raw_dihed_param_t> tori_list("charmm tori parameters");
  esmd::list<charmm_raw_dihed_param_t> impr_list("charmm impr parameters");
  esmd::list<charmm_raw_nonb_param_t> nonb_list("charmm nonb parameters");

  char buffer[INP_MAX_LINE_LEN];
  enum inp_sections cur_sec = PREAMBLE;
  while (fgets(buffer, INP_MAX_LINE_LEN, finp)) {
    char *line = inp_trim(buffer);
    //puts(line);
    if (line[0] == '!')
      continue;
    if (!strncmp(line, "BOND", 4)) {
      cur_sec = BOND;
    } else if (!strncmp(line, "ANGL", 4) || !strncmp(line, "THET", 4)) {
      cur_sec = ANGLE;
    } else if (!strncmp(line, "DIHE", 4) || !strncmp(line, "PHI", 3)) {
      cur_sec = TORI;
    } else if (!strncmp(line, "IMPR", 3) || !strncmp(line, "IMPH", 4)) {
      cur_sec = IMPR;
    } else if (!strcmp(line, "CMAP")) {
      cur_sec = CMAP;
    } else if (!strncmp(line, "NONB", 4)) {
      cur_sec = NONB;
    } else if (!strncmp(line, "HBOND", 5)) {
      cur_sec = HBOND;
    } else if (*line) {
      if (cur_sec == BOND) {
        charmm_raw_bond_param_t *bond = bond_list.get_slot();
        sscanf(line, "%s%s" REP2("%lf"), bond->itype, bond->jtype, &bond->kr, &bond->r0);
      } else if (cur_sec == ANGLE) {
        charmm_raw_angle_param_t *angle = angle_list.get_slot();
        int n = sscanf(line, "%s%s%s" REP4("%lf"), angle->itype, angle->jtype, angle->ktype, &angle->ktheta, &angle->theta0, &angle->kub, &angle->s0);
        if (n == 5)
          angle->kub = angle->s0 = 0;
      } else if (cur_sec == TORI) {
        charmm_raw_dihed_param_t *tori = tori_list.get_slot();
        int n = sscanf(line, "%s%s%s%s" REP3("%lf"), tori->itype, tori->jtype, tori->ktype, tori->ltype, &tori->kphi, &tori->n, &tori->phi0);
      } else if (cur_sec == IMPR) {
        charmm_raw_dihed_param_t *impr = impr_list.get_slot();
        int n = sscanf(line, "%s%s%s%s" REP3("%lf"), impr->itype, impr->jtype, impr->ktype, impr->ltype, &impr->kphi, &impr->n, &impr->phi0);
      } else if (cur_sec == NONB) {
        double ign, ign14;
        charmm_raw_nonb_param_t *nonb = nonb_list.get_slot();
        int n = sscanf(line, "%s" REP6("%lf"), nonb->itype, &ign, &nonb->emin, &nonb->rmin, &ign14, &nonb->emin14, &nonb->rmin14);
        if (n == 4) {
          nonb->emin14 = nonb->emin;
          nonb->rmin14 = nonb->rmin;
        }
      }
    }
  }
  ffpar->bond = bond_list.extract();
  ffpar->angle = angle_list.extract();
  ffpar->tori = tori_list.extract();
  ffpar->impr = impr_list.extract();
  ffpar->nonb = nonb_list.extract();

  ffpar->nbond = (int)bond_list.size();
  ffpar->nangle = (int)angle_list.size();
  ffpar->ntori = (int)tori_list.size();
  ffpar->nimpr = (int)impr_list.size();
  ffpar->nnonb = (int)nonb_list.size();
  return ffpar;
}

#ifdef TEST_CHARMM_INP
int main() {
  memory_init();
  FILE *f = fopen("../par_all27_prot_lipid.inp", "r");
  charmm_ffparam_t *ffpar = read_inp(f);
  ARR_FOR(v, ffpar->bond, ffpar->nbond) {
    printf("%s %s %f %f\n", v->itype, v->jtype, v->kr, v->r0);
  }
  ARR_FOR(v, ffpar->angle, ffpar->nangle) {
    printf("%s %s %s %f %f %f %f\n", v->itype, v->jtype, v->ktype, v->ktheta, v->theta0, v->kub, v->s0);
  }
  ARR_FOR(v, ffpar->tori, ffpar->ntori) {
    printf("%s %s %s %s %f %f %f\n", v->itype, v->jtype, v->ktype, v->ltype, v->kphi, v->n, v->phi0);
  }
  ARR_FOR(v, ffpar->impr, ffpar->nimpr) {
    printf("%s %s %s %s %f %f %f\n", v->itype, v->jtype, v->ktype, v->ltype, v->kphi, v->n, v->phi0);
  }
  ARR_FOR(v, ffpar->nonb, ffpar->nnonb) {
    printf("%s %f %f %f %f\n", v->itype, v->eps, v->sig, v->eps14, v->sig14);
  }
  memory_print();
}
#endif
